NASA Contractor Report 202314 
ICOMP-96-13 


O/^Z VJ8 


Simulation of Two-Fluid Flows by the 
Least-Squares Finite Element Method 
Using a Continuum Surface Tension Model 


JieWu 

Institute for Computational Mechanics in Propulsion 
Cleveland, Ohio 

Sheng-Tao Yu 
NYMA, Inc. 

Brook Park, Ohio 

Bo-nan Jiang 

Institute for Computational Mechanics in Propulsion 
Cleveland, Ohio 


December 1996 


Prepared for 

Lewis Research Center 

Under Cooperative Agreement NCC3^483 




National Aeronautics and 
Space Administration 



Simulation of Two-fluid Flows by the 
Least-Squares Finite Element Method Using a 
Continuum Surface Tension Model 


Jie Wu* Sheng-Tao Yu*and Bo-nan Jiang* 
NASA Lewis Research Center, Cleveland, OH 44135 


ABSTRACT 

In this paper a numerical procedure for simulating two-fluid flows is presented. This 
procedure is based on the Volume of Fluid (VOF) method proposed by Hirt and 
Nichols [1] and the continuum surface force (CSF) model developed by Brackbill, 
et al. [2]. In the VOF method fluids of different properties are identified through 
the use of a continuous field variable (color function). The color function assigns 
a unique constant (color) to each fluid. The interfaces between different fluids are 
distinct due to sharp gradients of the color function. The evolution of the interfaces 
is captured by solving the convective equation of the color function. The CSF 
model is used as a means to treat surface tension effect at the interfaces. Here a 
modified version of the CSF model, proposed by Jacqmin [3], is used to calculate 
the tension force. In the modified version, the force term is obtained by calculating 
the divergence of a stress tensor defined by the gradient of the color function. In its 
analytical form, this stress formulation is equivalent to the original CSF model [2]. 
Numerically, however, the use of the stress formulation has some advantages over 
the original CSF model, as it bypasses the difficulty in approximating the curvatures 
of the interfaces. 

The least-squares finite element method (LSFEM) [4] is used to discretize the gov- 
erning equation systems. The LSFEM has proven to be effective in solving incom- 
pressible Navier-Stokes equations and pure convection equations, making it an ideal 
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candidate for the present applications. The LSFEM handles all the equations in 
a u nifi ed manner without any additional special treatment such as upwinding or 
artificial dissipation. 

Various bench mark tests have been carried out for both two dimensional planar and 
axisymmetric flows, including a dam breaking, oscillating and stationary bubbles 
and a conical liquid sheet in a pressure swirl atomizer. 


1 Introduction 


Multi-fluid flows exist in many engineering problems. Examples of such flows include 
injection molding, metal casting, crystal growth and spray atomization, etc. At the 
interface of different fluids, surface tension exists as a result of uneven molecular 
forces. The interface behaves in a way similar to a thinly stretched membrane. 
The prediction of the evolution of the interface and the treatment of the interface 
conditions have been a challenging task for numerical simulations. 

The pressure jump across the interface is related to the surface tension coefficient 
a (only constant a is considered), the curvature of the interface k, and the viscous 
stress tensor r® by the Laplace’s formulas [5]: 

(P 1 -P 2 + = K.-fc - Tl ik )n k (1) 

where n is the unit normal of the interface, and the subscripts 1 and 2 denote the 
two different fluids. 

Most existing numerical methods for multiphase/ffee-surface flows fall into two cat- 
egories : (1) those which use a fixed grid; and (2) those which allow the grid to 
deform in time so that it remains surface-intrinsic. In the first category the compu- 
tational grid is fixed throughout the calculation. An additional variable is used to 
identify the interface (front). Examples of such methods axe the Marker and Cell 
(MAC) method proposed by Harlow and Welch [6] and the Volume of Fluid (VOF) 
method by Hirt and Nichols [1]. The MAC method uses massless marker particles 
which travel with the fluid to trace the fluids and the interface. The VOF method 
modifies the MAC method by replacing the discrete marker particles with a contin- 
uous field variable (color function). This function assigns a unique constant (color) 
to each fluid and has a sharp gradient at fluid interfaces. Numerical methods in this 
category are sometimes referred to as “front capturing” methods. Such methods 
possess great flexibility in handling large deformations and topological changes, as 
demonstrated by Daly [7] and Harlow and Shannon [8]. The most difficult task with 
the front capturing approach is to accurately identify the interface and to impose 
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the interface condition (1), such as exemplified by the elaborate work of Daly [9]. 
This difficulty can be alleviated by using a continuum surface force (CSF) model 
proposed by Brackbill, et al. [2]. When the CSF model is used, the interface condi- 
tion (1) is implied in the momentum equations rather than explicitly imposed, thus 
the location of the interface is no longer explicitly required in the calculation. The 
computer implementation of the CSF model is therefore relatively simple compared 
with other approaches. The combination of the VOF method and the CSF model 
has been used by a number of authors to simulate multi-phase phenomenon involv- 
ing surface tension and complex topological changes [10, 11]. A detailed description 
of the CSF model will be given in the next section. 

For methods in the second category [12, 13, 14, 15, 16], which are referred to as 
“front tr ackin g” methods, imposing the interface condition (1) is easy compared 
with the first, because the interface always coincides with mesh sides. However, it 
requires frequent updating of the computational mesh, which can be a complex and 
time-consuming procedure. In particular, it encounters severe difficulty when the 
flow experiences severe distortions and complex topological changes. 

Another approach which can be regarded as a combination of the above is the front- 
tracking method introduced by Unverdi and Tryggvason [17]. This approach uses 
a fixed, structured grid to represent the flow field. A separate, unstructured grid 
is used to represent the interface. The interface is explicitly tracked and kept at a 
constant thickness of the order of the mesh size. This ensures that the interface will 
not be smeared by numerical diffusion. Much success has been achieved in solving 
a variety of two-fluid flow problems using this approach [18, 19]. The difficulty with 
this approach is the handling of complex topological changes. 

In this paper numerical solutions to problems involving two immiscible fluids are 
sought. A large number of such problems deal with the interaction between a liquid 
and air, which are often simplified as free-surface problems. In the free-surface 
formulation the flow equations axe solved only for the liquid, and zero traction is 
assumed on the interface. In the present calculation such problems are treated 
as true two-fluid cases. The VOF approach and a modified version of the CSF 
model, proposed by Jacqmin [3], axe used in the simulation. The LSFEM [4] is 
used to discretize the governing system of equations. The LSFEM has proven to 
be effective in solving the Navier-Stokes equations for incompressible flows [4, 20] 
and pure convection equations [21, 22], m a kin g it an ideal candidate for the present 
applications. Some preliminary two dimensional calculations for two-fluid problems 
have been reported in [22]. The contents of this paper axe arranged as follows: In 
the next section the CSF model will be elaborated. In Section 3 the governing 
equations and the discretization procedure will be presented. Some two dimensional 
and axisymmetric numerical tests will be presented in Section 4. The first test 
case is a two dimensional dam breaking problem. The numerical results by the 
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present approach are compared with experimental data and with results obtained 
by other approaches. Good agreement is observed. The second test case deals with 
an oscillating bubble for both two dimensional and axisymmetric cases. The third 
case is a two-liquid interface problem in a closed box. The last test case simulates the 
formation of a liquid sheet in a pressure swirl atomizer (a Simplex nozzle). Finally, 
some conclusions are given in Section 5. 


2 The continuum surface force model 


The continuum surface force (CSF) model was introduced by Brackbill, et al. [2] 
as a means of treating the surface tension effect in the VOF method. By using 
this model, the difficulty in imposing the interface constraints (1) was alleviated. 
The basic idea of the CSF model is to regard the interface between two fluids as 
a transition region with a finite thickness, instead of a zero-thickness membrane. 
The surface tension effect is interpreted as a continuous body force spread across 
the transition region, which acts as a source term in the momentum equations. By 
using the CSF model, the interface condition (1) no longer needs to be explicitly 
imposed, as it is already implied in the momentum equations. The body force in 
the CSF model is obtained through some differential operations on the spatially 
continuous color function, as follows: 

,= r vc 

where f is the body force, a and « are the surface tension coefficient and the cur- 
vature of the interfaces, the same as in (1), C is the color function, and [C] denotes 
the jump of C across the interface. The curvature k in the above formulation is 
calculated from: 

/c=-(V-n) (3) 

where n is the unit normal to the surface, which is obtained from: 

VC 

n “|vc| 



( 2 ) 


Jacqmin [3] derived the CSF model through the analysis of tension energy. He 
pointed out that the surface tension force can be expressed as the divergence of a 
stress tensor which is uniquely defined by the gradient of the color function. Let 
Tij denote the the stress components in a Cartesian coordinate system. They are 
related to the color function C as : 

<r 1 fdCdC dC dC\ 

Tii ~ [c]|vc| ” bij dx k lhr k ) 
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(5) 


<J 


[C\ 


( 1 dC dC 
\JVC| fo,- dxj 


- 6a |VC| 


where 8ij is the Kronecker delta. 

The volumetric body force caused by the surface tension is expressed as: 

f _ 

Jt d Xj 


(6) 


It can be easily verified that f defined by Eqs. (5) and (6) is analytically equivalent 
to that of the original CSF model given by Eqs. (2)-(4). The advantage of using the 
formulation given by Eqs. (5) and (6) instead of the original CSF model is that Eqs. 
(5) and (6) do not require the explicit calculation of the normalized gradient term in 
the right hand side of Eq. (4), whose definition is not cleax when the denominator 
| VC| approaches zero. In Eqs. (5) and (6) both r and f are well defined in the whole 
domain; and naturally vanish when |VC| becomes zero. An additional advantage 
of using the stress tensor t is that it can be regarded as part of the momentum 
flux. In many numerical procedures t can be used directly and there is no need to 
calculate f . 

It is worth pointing out that the stress tensor r also describes more directly the 
fluid physics at the interface than the force term f. To illustrate this argument, let 
1 and 2 denote, respectively, the directions along and normal to the interface in two 
dimensions (Figure 1). We have: 

(VC)! = 0, (VC) 2 = |VC| 


The stress components are 
ure 1): 

Tn = 


(for clarity the definition of stress is also given in Fig- 



t 22 — r i2 — r 2i — 0 


From the above we have: 



= cr 


This indicates that when the interface thickness approaches zero, the stress tensor 
defined by Eq. (5) also approaches the true surface tension. The fact that Tn is 
always positive in the transition region reflects that the interface is like a stretched 
membrane, t is also useful in situations when the surface tension coefficient a is no 
longer constant. 
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Figure 1: The continuum surface and the stress components 

It should be pointed out that Lafaurie, et al. [11] introduced a capillary pressure 
tensor T, and used it as a correction to the momentum stress tensor. T is defined 
as: 

T = a(I-n<g>n)|VC| (7) 

where n is the unit normal to the surface as defined by Eq. (4). It is easy to verified 
that t in Eqs. (5) and (6) and T in Eq. (7) are also analytically equivalent except 
for a minus sign. Again the advantage of using t is that it bypasses the calculation 
of n. 


3 The Governing Equations and the Discretiza- 
tion Procedure 


The governing equations for the two-fluid problems are the standard Navier- Stokes 
equations. Incompressibility is assumed for both fluids. In order to use the Co 
element in the LSFEM formulation, the governing equations need to be re-written 
as a first-order system. Here we choose the following velocity- pressure- vorticity 
form: 

V • u = 0 (8) 

+ p (u • V) u + Vp + /i ( V x u>) 

- (V/x • Vu + Vp- (Vu) T ) = f (9) 

in which p is the density, u is the velocity vector, p is the pressure, u> is the vorticity 
vector, p is the dynamic viscosity, and f is the body force, which generally consists 
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of the surface tension effect given by Eq. (3) and the gravitational force. The terms 
in the last pair of brackets on the left hand side represents the effect of non-uniform 
viscosity. 

The vorticity u> is defined by the velocity u as: 

OJ = V x u (10) 

Prom the above we immediately have: 

V • u> = 0 (11) 

The above equation system, Eqs. (8)-(ll), has been used by Jiang et al. [23] in the 
theoretical analysis of Navier-Stokes equations systems; and has been extensively 
used as the basis of LSFEM calculations [4, 20, 24]. 

The fluids are identified by the different value of the color function C, which is 
convected by the flow field: 

+ (u • V)C = 0 (12) 

Ol 

Some numerical examples in this paper deal with axisymmetric (non-swirling and 
swirling) flows, which require the above equations to be recast into a cylindrical 
coordinate system. These equations will be given in Appendix A. 

Fluid properties such as the density and the viscosity are assumed to be distributed 


in the same manner as C, i.e., 

P = + ( 13 ) 

C/2 — C/i 

p = * + < c - Ci > (14) 

The governing equations, Eqs. (8)-(12), are first discretized in time. The backward 
Euler difference scheme is used: 

V ■ u n+1 = 0 (15) 

pn+lU " + A^~ + ^ ( U ” +1 ' V ) Un+1 + Vp ” +1 + pn+1 ( V X Wn+1 ) 

- • Vu n+1 + V/x n+1 • (Vu" +1 ) T j = r +1 (16) 

U) n+1 = V x u n+1 (17) 

- + (u n+1 • V)C n+1 = 0 (18) 
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The superscript ‘n’ denotes the previous time step and ‘n + 1’ denotes the current 
time step. 


The above equations, Eqs. (15)-(18), are further linearized. To ensure time-accuracy 
for time-dependent problems linearization iterations are performed within each time 
step (typically three iterations are used). For steady state problems only one lin- 
earization iteration is performed in each time step. Let r denote the number of 
linerization iterations, the linearization is carried out as follows: 


V • (u n+l ) = 0 
V )\r] 

( «+i\ ^ + (p'+i) (7u” +1 ) v')(u" +1 ) 

+ V (p” +, ) W + ('‘* + %- 1| V > x (“* +1 ) m 
= ( r+1 )|'-i| + (' , " +l ) l '-, l li 

+ ( V (' t ” +1 ),.-, l • V ( u * +1 )['- 1| + v • ( v K 1 ),'.,,)' 

(C n+1 ) 

aT 


(19) 


+ ((“" +l ) w - v )( cm+1 ) l '] = ^ 


zii 

At 


(20) 

(21) 

( 22 ) 


In the above (C” +1 ) (0] = C n , (^ n+1 ) [0l = H n and (u n+1 ) [0] = u n 


[o] 


no] 


At this stage the standard LSFEM procedure [4] is introduced for spatial discretiza- 
tion. We note that the solution of the flow field ((u n+1 )j r ], (p w+1 )[ r ], (u> n+1 )[ r ]) does 
not require the knowledge of (C' n+1 )[ T ], Thus for each linearization iteration, Eqs. 
(19-21) are first solved. Then the newly obtained (u n+1 )[ r j is used in Eq. (22) to 
solve for (C n+1 )[ r ]. We use reduced integration (one point Gauss quadrature for 
bi-linear elements) in the solution of Eqs. (19)-(21). The reason for so doing is that, 
as pointed out by Jiang and Povinelli [4], the LSFEM with Gaussian quadrature is 
equivalent to the collocation least-squares method; and the use of reduced integra- 
tion makes the total number of collocation points compatible with the total number 
of unknowns and has proved to result in better accuracy [21]. For Eq. (22), however, 
we have used exact numerical integration (2x2 Gauss points for bi-linear elements), 
because using reduced integration sometimes produces oscillatory solutions for the 
color function. Such oscillations will further result in unacceptable errors in the 
tension force, which is obtained by twice differentiating the color functions. In order 
to avoid excessive numerical diffusion and to maintain time accuracy, the time step 
is chosen so that the corresponding CFL number is less than or close to unity for 
true transient problems. If only steady state results are of interest, much larger 
time steps are used. A Jacobi-preconditioned conjugate gradient (JCG) method has 
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been used to solve the algebraic equations. Since the LSFEM always results in a 
symmetric, positive-definite matrix, the JCG method is guaranteed to converge to 
the correct solution. Since the JCG iterations can be carried out in an element- 
by-element fashion such that no global matrix needs to be stored [25], very large 
problems can be solved with moderate computer memory requirements [20]. 

The body force term (f n+1 ) [r _ 1] in Eq. (20) is calculated using Eqs. (5) and (6). The 
finite element solution of r obtained from Eq. (5) is discontinuous across element 
boundaries. In order to use Eq. (6) to obtain the body force f, the r field needs 
to be spatially continuous. Such a continuous field can be obtained by treating 
the dis continuous field with various recovery procedures (see [26] and the references 
therein). In this paper only bi-linear finite element is used for the discretization 
and a simple averaging method [27] is used for the recovery procedure. For higher 
order elements, the superconvergent patch recovery (SPR) [26] should be used. The 
recovery procedure is also used to obtain a continuous f field. 


4 Numerical Examples 

4.1 A Broken Dam Problem 

This problem has been used by many as a test case for simulating free-surface 
problems. Experimental data for this case are available [28]. Here the problem is 
solved as a two-fluid problem involving both the water and air. Zero surface tension 
and slippery walls (left and bottom sides) are assumed. On the top and right sides 
zero pressure is imposed. In addition, if inflow is detected on the top and right 
sides, the density of air is imposed. The computational domain is 2 units high and 
6 units long. Ini tially, water occupies a 1 x 1 area at the bottom left corner. The 
computational mesh consists of 120 x 40 uniform bilinear quadrilateral elements. 
The time step is A t = 0.05. The nondimensionalized gravitational acceleration, g, 
is taken to be unity. The viscosity is set at 3.05 x 10“ 5 for water (which is the same 
choice as in [29]) and 3.05 x 10 -8 for air. The densities for water and air are 1 and 
0.001, respectively. 

Figure 2 shows free surface profiles (a contour line of density at p = 0.5) and the 
pressure contours at various times. The calculated water front location and water 
column height are compared with the experimental data in Figure 3. It can be seen 
that the calculated results are in good agreement with the experimental data [28] 
and with the calculation results in [29]. 
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Figure 2: 
bottom: t 


A broken dam, free surface profile and pressure contours. From top to 
= 0.5, 1.0, 1.5, 2.0, 2.5, 3.0. 
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Figure 3: A broken dam, front location and water column height, solid line — 
calculated results; dots — experimental data. 

4.2 An Oscillating Bubble 

In this example the oscillation of a bubble due to surface tension is studied. Both 
two-dimensional planar and axisymmetric cases are considered. The two dimensional 
case was simulated by Fyfe et al. [13] using a Lagrangian approach. The compu- 
tational domain in this calculation is 4 x 12 units with the symmetry plane/axis 
along the longer side. The computational mesh consists of 40 x 120 uniform bilin- 
ear quadrilateral elements. Initially the shape of the bubble is elliptical, given by 
x 2 /4 -f z 2 = 1. The density of the fluid inside and outside the bubble is 1.5 and 0.5, 
respectively. The dynamic viscosity is 0.01 for both fluids. A time step of At = 0.05 
is used in both cases. After the calculation reaches f = 15.0, viscosity is increased to 
1.0 to allow the solution to converge to steady state. When steady state is reached 
the sizes of the bubbles are 1.40 and 1.58, respectively, while the theoretical values 
are 1.41 and 1.59. At steady state, it is verified that the Laplace’s formulas Eq. (1) 
for the pressure jump is satisfied. 

In Figs. 4 and 5 the interface shapes at different times are shown for two dimen- 
sional and axisymmetric cases, respectively. The corresponding time history of the 
interface locations on the x -axis is shown in Figure 6. The oscillation periods are 
estimated (based on the first two cycles) to be 7.6 for the two-dimensional case and 

7.2 for the axisymmetric case. Theoretical solutions of the oscillation periods are 
available from a linear analysis for small amplitudes (see [13, 19] for details and 
references), which for this example axe 6.1 for the two-dimensional case and 6.3 for 
the axisymmetry case. Since in the calculations both nonlinear and viscous effects 
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Figure 4: An oscillating bubble, interface profile for two dimensional case, 
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 40, respectively 
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Figure 5: An oscillating bubble, interface profile for axisymmetric case, t 

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 40, respectively 
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Figure 6: An oscillating bubble, time history of the interface location on the x -axis, 
top: two dimensional case; bottom: axisymmetric case. 
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are present, the theoretical solutions can be used only as references. The current 
numerical simulation predicts higher values of the oscillation periods than Unear 
theory, which is consistent with the experience of Fyfe et al. [13]. 


4.3 Two-liquid Interface Problem 


This example is due to Tezduyar, et al. [16]. In this problem, a two dimensional 
0.8 x 0.6 box is fully occupied by two liquids of equal volume. The density of the 
fluids are 1.0 and 2.0, respectively. The gravitational acceleration is 0.294. The 
Ughter Uquid is placed on top of the heavier one. Initially the interface is a straight 
line, with a slope of 0.25. The dynamical viscosity is 0.001 (the same for both 
Hquids). Zero surface tension is assumed. The computational mesh consists of 
21 x 42 uniform bilinear elements. A time step of At = 0.25 is used throughout 
the calculation (which corresponds to a CFL number close to 1 when the velocity 
reaches its maximum at approximately t = 6). SUppery wall conditions are used 
on all four sides. Figure 7 shows the interface profile and velocity field at different 
times. Figure 8 shows the time history of interface location on the left and right 
hand side walls (relative to and normalized by the average height). The results agree 
well with those presented in [16]. 


4.4 Flow in a Pressure Swirl Atomizer 

T his example deals with the swirUng flow in a pressure swirl atomizer (a Simplex 
nozzle). Simplex nozzles are widely used as fuel atomizers in gas turbine engines. 
The formation of a conical liquid sheet is characteristic of such nozzles. Prediction 
of the thickness of the liquid sheet and the exit angle for various inlet conditions < 
provides helpful information in the design process. The computational domain for 
the current simulation is shown in Figure 9. The liquid enters to the swirl chamber 
though a number of inlet slots. For the present calculation, axisymmetry is assumed. 
Gravitational effects are ignored. The computational mesh, shown in Figure 10, 
consists of 3600 bilinear elements with 3741 nodes. Initially the whole domain is 
filled with motionless air. At t > 0 liquid begins to stream into the domain through 
the inlet slots. At the inlet, values of the liquid density and velocity are prescribed. 
The prescribed values are: u r = -4.8, u $ = 8.0, u z = 3.6, p = 0.1 The density of 
the air is p = 0.01. The dynamic viscosity for both fluids axe fi = 0.005. For nodes 
on the free boundary, the density of air is imposed if inflow is detected, nothing 
is imposed if outflow is detected. Slippery boundary conditions are imposed on 
portions of the wall; and non-slip conditions are imposed elsewhere on the wall, as 
shown in Figure 9. The pressure is specified on the top left comer of the domain. A 
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Figure 7: Two- liquid interface probli 
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 


, interface profile and velocity field, t = 
, 16, 17, 18, 19, respectively 
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Dimensionless wave height 



Figure 8: Two-liquid interface problem, time history of interface location on the left 
and right hand side walls 
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constant time step At = 0.1 is used throughout the calculation. As only the steady 
state solution is of interest, a linearization iteration is performed only once in each 
time step. Figure 10 shows the steady state interface profile. The formation of the 
conical liquid sheet is successfully predicted in this simulation. This is a preliminary 
simulation and no attempt has been made to compare computational results with 
measured data. At the moment, no experimental data are available for comparison. 



Figure 9: Flow in a pressure swirl atomizer, the problem definition 
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Figure 10: Flow in a pressure swirl atomizer, the computational mesh and the 
interface profile 

5 Concluding Remarks 


In this paper a numerical procedure based on a CSF model and the LSFEM is pre- 
sented for two-fluid flow problems. Numerical tests carried out on a number of two 
dimensional planar and axisymmetric flows indicate that this approach is capable 
of simulating such flow phenomena. The present approach has the advantage in the 
ability to handle complex topological changes such as the breakup of a liquid jet. 
Another advantage is the ability to handle complex geometrical configurations in 
practical engineering environments, since the LSFEM is based on totally unstruc- 
tured grids. Currently the interface is typically spread over 3-5 grids. A change 
of the interface thickness is also observed in the calculations. The analysis by Haj- 
Hariri [30] indicates that the migration velocity is strongly affected by the smearing 
of interface. Thus, to ensure good accuracy of the simulation a fine grid would have 
to be used. Improvement of the accuracy of the numerical scheme is expected by 
introducing higher order discretization schemes in both time and space. Our tests 
[22] indicate that the discontinuity can be modeled within 3 grids, even when the 
interface undergoes very large deformation due to a highly vortical flow field. 
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A The Governing Equations in Cylindrical Co- 
ordinate System 

The Navier-Stokes equations given by Eqs. (8)-(9) can be written in a cylindrical 
coordinate system as: 
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The definition of u>, Eq. (10), is given as: 
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The compatibility condition, Eq. (11), in this system becomes: 
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The stress components due to surface tension are: 
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in which |VC| = • 

The body force components due to the above stress tensor is: 
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(39) 

(40) 


When the flow is axisymmetric all the ^-derivatives become zero. Further if the flow 
is non-swirling, uq becomes zero. For such cases the governing equations are further 
simplified. 
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